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Hadronic matrix elenrents of proton decay are essential ingredients to bridge the grand unification 
theory to low energy observables like proton lifetime. In this paper we non-perturbatively calculate 
the matrix elements, relevant for the process of a nucleon decaying into a pseudoscalar meson and 



O ; an anti-lepton through generic baryon number violating four-fermi operators. Lattice QCD with 

,—, . 2+1 flavor dynamical domain- wall fermions with the direct method, which is direct measurement 

of matrix element from three-point function without chiral perturbation theory, are used for this 
tH- " study to have good control over the lattice discretization error, operator renormalization, and chiral 

tj- ■ extrapolation. The relevant form factors for possible transition process from an initial proton or 

t>; 

_il . neutron to a final pion or kaon induced by all types of three quark operators are obtained through 

~f~, \ three-point functions of (nucleon)- (three-quark operator)- (meson) with physical kinematics. In this 

study all the relevant systematic uncertainties of the form factors are taken into account for the 
K> . first time, and the total error is found to be the range 30%-40% for n and 20%-40% for K final 

states. 
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I. INTRODUCTION 

Proton decay is a smoking gun evidence of physics 
a natural outcome of Grand Unified Theories (GUTs) 
baryon number changing interactions mediated by the heavy new particles. Dominant modes 
are of X and Y gauge boson exchange for non supersymmetric (SUSY) GUTs and of color- 



aeyond the standard model and is 
l|, |2j. The process occurs through 



triplet Higgs multiplet for SUSY GUTs 



flEl 



4J. Recent SuperKamiokande experiments report 



the bound on proton partial lifetime, for instance, r > 8.2 x 10 33 year for the p — > e + 7r° 
channel [5|, |6j, which is typical forgauge boson exchange, or r > 2.3 x 10 33 for p — > K + u [7J 
and r > 1.6 x 10 33 for p — > K°£i + 8[, both of which are favored for some SUSY GUTs. There 
have been many arguments of a constraint on proton lifetime from various types of GUT 
model so far (see a comprehensive review [9j and reference therein). In order to constrain 
the parameter space in GUT models with a reliable bound, a removal of every theoretical 
uncertainty is highly desirable. One of the important elements, which can be made less 
uncertain from the current knowledge, is the hadronic contribution to proton decay matrix 
elements. Lattice QCD calculation can lead to reducing the uncertainties in hadronic matrix 
element of a nucleon decaying into a pseudoscalar meson, and thus it can provide relevant 



information for the proton lifetime bound and help experimental plans for the future [10] . 

The estimate of proton decay matrix elements in lattice QCD has been significantly 
improved by removing systematic errors, one by one, since the first attempts in 1980s 
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131. A decade ago using Wilson fermion action with quenched QCD, JLQCD collaboration 
141 ] performed an extensive calculation of proton decay matrix elements with both the 
"direct 1 method, which is a direct measurement of matrix element from three-point functions, 
and the "indirect 1 method, which is an effective estimate through low-energy constants in 
tree-level chiral perturbation theory, calculated with two-point functions. Few years later 
JLQCD and CP-PACS joint collaboration carried out a continuum extrapolation of the 
low energy constants for the indirect analysis [l5| to remove the uncertainty due to large 
discretization error. Within the "direct" method, RBC collaboration 16[ performed the 
analysis using domain-wall fermions (DWFs) and a non-perturbative renormalization, where 
thanks to almost exact chiral symmetry of the DWFs the discretization error of 0(a) is 
essentially removed and the error of renormalization factor associated with the use of lattice 
perturbation theory was also eliminated. The RBC collaboration and later RBC/UKQCD 



collaborations extended the DWF calculation of the "indirect" method, in which they used 



the quenched approximation [16| as well as unquenched one containing active two (u and 
d) and three (adding s) dynamical light flavors [17J. In this way, one of the uncontrolled 
systematic uncertainty coming from quenched approximation was removed. 

A striking, but perhaps not surprising outcome of the comparison of the results from 
direct and indirect calculations, though performed only with quenched approximation so far, 
is that the indirect method could overestimate the matrix elements by a factor of about 

n 

two [16]. Therefore the last step is to perform the direct calculation with the three flavor 
DWF simulation with non-perturbative renormalization of the operators, to remove all the 
uncontrolled systematic uncertainties in the existing calculations. 

In this paper we provide the non-perturbative value of proton decay matrix elements 
using the direct method with the dynamical, Nf = 2 + 1 (degenerate u, d and physical s 
quarks) flavor lattice QCD with DWFs. The DWF ensemble for Nf = 2 + 1 at the lattice 



cutoff a" 1 ~ 1.7 GeV with 300-700 MeV pion masses [18| in RBC/UKQCD collaboration 
are used for this purpose, and thus this enables us to evaluate hadronic matrix elements 
including almost all systematic errors on the lattice. 

This paper is organized as follows. In section [III we explain the definition and property 
of the matrix elements as well as their relation to the proton partial decay width. The 
method to extract the matrix elements from three-point function on the lattice is expressed 
in section IHIl and in section IIVI we present our setup and the detailed analysis to obtain 
the matrix elements and evaluate their systematic uncertainties. Section [V] is devoted to 
summary and outlook. 

II. PROTON DECAY MATRIX ELEMENT 

A. Effective Lagrangian and matrix element 

Baryon number violating operators appearing in the leading low-energy effective Hamilto- 
nian are constructed by possible combination of dimension-six (three quarks and one lepton) 
operators to be SU(3) color singlets and SUl(2)x Uy(l) invariant. Following the notation 



of |l9H2l|. four-fermi operators are expressed as 

{ l d = (DlUiUq^l^e^, (1) 

Oil = {CA P )MM) R e ijk e a ^ (2) 

oitL = (C, A(gr,^V^, (3) 

OIL = {UlD{) R {Ull d ) R e^\ (4) 

with generic lepton field Z, and quark field of left-handed part q and right-handed part U 
and D as up and down type. The indices a, b, c, d denote the generation number of fermion, 
i,j, k denote color SU(3) indices, and a,/3,j,5 are SU(2) indices. The inner product is de- 
fined as (x, h)r/l = ^ t CPrilV which has charge conjugation matrix C and chiral projection 
Pr/l- The baryon number violation (but preserving B — L number) in GUT models is gener- 
ally expressed as low-energy effective Hamiltonian with the above six-dimension operators. 
Leading term of effective Hamiltonian at low energies reads 

^=J2C I [(qq)(ql)] I ^) + --- = - £ ^[TOjV) + ■ • • , (5) 

I I 

where C 1 = C \fi) is the Wilson coefficient with renormalization scale \x for each operator 
i corresponding to Eq.flTJ-fllD. The details of the (SUSY) GUT is all captured in the coef- 
ficients C 1 ^). Ellipsis means the higher order operators which are suppressed by inverse 
power of heavy mass scale. The index % within the three-quark operator O qqq and lepton I 
in Eq.(J3]) distinguishes the type of fermion and lepton (flavor and chirality). Three-quark 
operator is defined as 

OZ = (qq)rqr = e^(q' T CP r q^P ri q k , (6) 

with light quark flavor q = (u, d, s), where the color singlet contraction is taken. Dirac spinor 
indices are omitted in the above equation. In the following we may use simple notations 
for the three-quark operators as O vv . V and r" denote the chirality, either R or L and the 
bracket means the contractions among Dirac spinors. 

We calculate the transition matrix elements of the dimension-six operators with an initial 
nucleon (proton or neutron, N = p,n) state and a final state containing a pseudoscalar meson 
(P = (jr,K,rj)) and an anti- lepton (I) 

(P(p),l(q,sWO rr ']\N(k)) = vf(q,s)(P(<p)\O rr, \N(k,s)), (7) 



including three-dimensional momenta, p for final pseudoscalar, k for initial nucleon and 
q = p — k for final lepton which is determined from momentum conservation. Neglecting 
the electroweak interaction of the lepton, the amplitude (l(q, s)\l c \0) = vf(q, s) of the lepton 
part can be captured in the wave function of on-shell lepton state at momentum g*for spin s 
component. Furthermore the matrix element in the quark sector can be divided into relevant 
form factor Wo(q 2 ) and irrelevant one Wi(q 2 ) as 



tTT 



rr' 



(P(p)\O lL \N(k,s)) = P r W^ (q 2 )-^W± L tf) u N (k,s). 



m N 



(8) 



where Wq and W\ are defined for each matrix element with the three-quark operator renor- 



malized in MS NDR at scale /i, and are functions of square of four momentum transfer 
q = k — p. Using on-shell condition, the total matrix element as shown in Eq.(JTI) is given by 

id 



v?(q,s)(P(p)\0" \N(k,s)) = vt(q,s)P T , W^ (q 



rVF' i 



m N 



W[ l (q 2 ) u N (k,s) 



rr'. 



vf(q, s)Pr>u N (k, s)W^ 1 (0) + 0( mi /m N 



(9) 



with ifjivi = rriiVi and W\ ~ Wq 16|. Since — q 2 = m 2 is much smaller than nucleon mass 



squared in the case of / = e, z/, we set q 2 = and ignore the second term in Eq.([H]). Taking 
only the relevant form factor will be a good approximation even for / = /i, as vn^jm^ ~ 10% 
is smaller than the total error of Wq in this study. 

Once we obtain the relevant form factor Wq in lattice QCD, the partial decay width of 
the decay iV — > P + / is given by 

with the perturbative estimate of Wilson coefficient C 1 in the GUT models [9|. Note that 
renormalization scale dependence of Ci and Wq cancels out in their multiplication. 

We can reduce the number of types of matrix element by using Parity symmetry and 
isospin symmetry. The different chirality combinations of the matrix elements are related 



through the Parity transformation as 

(P;p\0 RL \N; I s) = l0 (P; -p\0 LR \N; -k, s), 
(P;p\0 LL \N; k, s) = l0 (P; -p\0 RR \N; -k, s), 



(11) 
(12) 



which indicates that four chirality combinations (rr') = (RL), (LL), (RL), (RR) are reduced 
to two different combinations, (rr') = (RL), {LL). In this paper V is fixed in a left-handed 



chirality, and thus we define Wq{ = Wq X . Under exchange-symmetry between u and d there 
is the relation between their isospin partners n and p, 



(n \(ud) r u L 
(7i + \(ud) r d L 
(K°\(us) r u L 
(K + \(us) r d L 
(K + \(ud) r s L 
(K + \(ds) r u L 
(r)\(ud) T u L 



\p) = (7r°\(du) r d L \n), 

\p) = -{7r~\(du) r u L \n), 

\p) = -(K + \(ds) r d L \n), 

\p) = -(K \(ds) r u L \n), 

\p) = -(K \(du) r s L \n), 

\p) = -(K \(us) r d L \n), 

\p) = -(v\(du)rd L \n), 



(13) 
(14) 
(15) 
(16) 
(17) 
(18) 
(19) 



involving negative sign coming from u and d exchange of interpolation operator of proton 
and of neutral pion. Furthermore in the SU(2) isospin limit there is an additional relation 
between Eq.flU and Eq.(fT4l): 



{n°\(ud) r u L \p) = V2(n + \(ud) r d L \p). 



Therefore there are twelve principal matrix elements we calculate in this paper. 



(20) 



III. CALCULATION SCHEME FOR THE FORM FACTORS 

To obtain the matrix element we make use of the ratio of three-point function of (nucleon)- 
(O rL )- (meson) and two-point function of nucleon and meson. Such a ratio is represented 
as 



Ra(t,t u to;p,P) 

E^^^^tr^^lJf^i^QO^^^Jf^to)^)] 

= E* A e<«ft-*><0| J? (aMO J? \x,t)\0) E,-tr[P 4 (0| Jf (x,t)Jf (0,t )|0)] 



(21) 



with interpolating field for pseudoscalar Jp and nucleon J p . Here we introduce the spatial 
momentum p = 2irn/ 'L a , n is integer vector taking from to L a — 1, L a is the spatial 
extension of the lattice. Superposition "gs" of interpolating operator denotes the type 
of source/sink function is gauge-invariant Gaussian smearing 22j. "tr" represents trace 
over spinor indices, and V is a spin projection matrix. In order to obtain the three-point 
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function in numerator including three-quark operator defined in Eq.flfJ]), we construct quark 
propagator with sequential source method. 

In this calculation we use n = (1,0,0), (1, 1,0) for meson momentum and zero nucleon 
momentum in three-point function. Zp >p indicates the amplitude of overlap of the interpo- 
lating field to on-shell state, 



(P(p)\Jf t(o)|0) = ^Zf(p), (22) 

(0|Jf(0)b(0,s)> = yfzfu p {m N ,s), (23) 

with Dirac spinor normalized by u p (k, s)u p (k, s') = 2mjq8 ss i being the momentum of the 
nucleon k = (0, 0, O^m^) and meson p = (p,iEp). Note that the operator of nucleon 
interpolating field is not uniquely determined, and actually we see that two possible proton 
operators formed as 

J p = e ijk {u iT C^d j )u k , e ljk {u iT C l4l ^)u k . (24) 



Since there is equivalence between proton and neutron matrix element as shown in Eq. (TT3|) - 
ffT9|) . we focus on proton case. Numerical comparison between the above two types of nucleon 
interpolating operator will be shown in the next section. 

Taking an enough large time separation of meson-operator (ti — t) and operator-nucleon 
(t — 1 ) in order to suppress the other excited state contamination than an asymptotic state 
of pseudoscalar and nucleon, the ratio leads to 

lim R 3 (t,tiMP,V) = R7 m (P^) = tr[m«(g 2 ) - -£- W?(q 2 ))}, (25) 

,t— to— too L TflAT J 



£]_— t,t— to 



m N 



where q 2 is the squared momentum transfer from the initial proton to the final pseudoscalar 
meson state q 2 = (k — p) 2 . We employ two different projection matrices V = P 4 or iP^jj 
with P 4 = (1 + 74) /2 by which we subtract parity-partner contamination (excited state) to 
extract form factor from three point function by solving linear equation of 

it>r yn >, P 4 ) = W T Q {q 2 ) - ^<(g 2 ), (26) 

m N 

RT ym (P^P^) = ^W[(q 2 ), (27) 

and thus Wq and W\ can be simultaneously obtained. Hereafter we concentrate on the 
relevant form factor Wq. 



IV. NUMERICAL CALCULATION OF THE PROTON DECAY FORM FACTORS 

A. Lattice setup 

We use gauge configurations with dynamical domain-wall fermions of 2 + 1 flavor at 
wasaki gauge action in lattice size 24 3 x 64 at = 2.13 which corresponds to a" 1 = 1.73(3) 



181 ] . This is the same ensemble as previous indirect method study [l7j . Boundary condition 
is periodic for gauge field, and spatially periodic and temporally anti-periodic for fermion 
field. We use four different unitary u, d quark masses in chiral extrapolation, and one 
unitary and one partially quenched strange-quark mass in the study of strange quark mass 
dependence for final K°' + kaon state. Explicit chiral symmetry breaking due to finite fifth 



181 ] (residual mass is m rcs ~ 3 x 1CT 



dimensional lattice as L s = 16 is small in this ensemble 
and indeed we take into account residual mass-shift in the chiral extrapolation. For later 
convenience let us introduce the quark mass fh which includes the additive renormalization 
due to the inexact chiral symmetry of the domain-wall fermions at a finite extension of the 
fifth dimension. We define 

fh = m + m res , (28) 

as the multiplicatively renormalizable mass with m in the lattice action, where residual 

n 

mass m res for the lattice used in this study has been calculated as m res = 0.003152(43) [18j. 
For the nucleon to pion matrix elements we have two controlling variables: rh u d for the 
degenerate u and d quark mass and the squared momentum transfer q 2 . For the nucleon to 
kaon matrix elements the strange quark mass fh s is an added controlling variable. 

When computing two-point and three-point function on the lattice, we employ parameters 
of gauge-invariant Gaussian smearing with its parameter (no-, cr) = (40, 5.0) for baryon 



source/sink and (no, cr) = (16, 3.0) for meson sink including AP 



smeared link gauge in 



which we use the parameter (N,c) = (12,0.4) as defined in [23j. Nucleon source-point 
to and pseudoscalar sink-point t\ for three-point function are set to be fixed as to = 5 
or 37 and t\ = 27 or 59, and operator in three-point function constructed by sequential 
source for meson sink moves its position t between to and ti with spacial momentum p = 
(7r/12,0, 0), (7r/12,7r/12,0). The detailed number of gauge ensembles and source positions 
employed here are shown in table |T] For measurement at m s ^ = 0.005 ensemble we add 
the statistics in different source and sink points in the same configurations separated in 40 



HMC trajectory skips. At m s ^ = 0.02, 0.03 the skipped HMC trajectory in the same source 
and sink point is 40 and in every source-sink combinations we separated 20 HMC trajectory 
in order to reduce the autocorrelation effect. At m s ^ = 0.01 the skipped HMC trajectory 
is 80 and separated trajectory for every source and sink combinations is 40. Note that at 
the lightest quark mass (m^f = 0.005) to increase statistics we further take average over 
two different source time slices, to = 5 and 37, for each gauge configurations (therefore the 
number of measurements is double of configurations) . 

The multiplicative renormalization factors to convert the lattice three-quark operators 
in Eq.( ll3p -( 1T9"]) into those in MS NDR scheme has been calculated through the RI/MOM 
non-perturbative renormalization [17| as 

£/(/i = 2GeV) LL = 0.662(10)(53), (29) 

U((i = 2GeV) RL = 0.665(8)(53). (30) 

The first error is statistical one and the second is systematic one (8% systematic error is 
estimated in [17J] as a truncation effect of perturbative expansion). 

In Figured] we show the effective mass plot of nucleon, pion and kaon two-point function 
which consists of denominator of Eq.( l2~lj) . This is constructed with data at t and t + 1, and 
we can observe the plateau region whose starting point is t — 5 for nucleon and t = 6 for 
pseudoscalar. It turns out that the asymptotic state of nucleon and meson state of three 
point function in Eq.( l2ip will be reached in the region t — 1 > 5 and t\ — t < 6 for temporal 
position of operator t. 

B. Measurement of form factor and kinematics 

Figures [2] and |3] show the form factor Wq of the p — > it channel obtained by using the 
asymptotic structure Eqs. ( 1261) and (1271) as a function of the time position t of the three-quark 
operator for each quark masses and momenta. The open and filled symbols correspond to 
results in two different nucleon interpolating operators, {q T C , y5q)q and {q T C , ~f4'-f^q)q respec- 
tively. To obtain the value of Wq, a simultaneous fit of these two effective Wq is performed 
at the plateau in the range 13 < t < 20, where the two Wq appear to be consistent and the 
contamination from the excited states dies out. The same range is used for all the parame- 
ters and all the matrix elements. Figures |4] and [5] show W for each channel as a function 
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TABLE I: Lattice parameters used in this paper. Fitting range of pion and Kaon mass is 6 < t < 23, 
and nucleon mass is 5 < t < 13. Lattice spacing is a" 1 = 1.733 GeV, which is same value as la] . 
Two columns of — q 2 are values of squared momentum transfer for each quark masses and two 
different spatial momentum squared, |p| 2 = (tt/12) 2 , 2(tt/12) 2 respectively. 



(m^jinf 3 ) (m^,mj al ) 



m, 



m K 



m N 



Y(GeV 2 



# configs. # meas. 



(0.005,0.04) 


(0.005,N/A) 


0.1897(5) 


N/A 


0.656(16) 


-0.129 


0.241 




(0.005,0.0343) 


N/A 


0.3131(5) 


N/A 


0.017 


0.325 




(0.005,0.04) 


N/A 


0.3322(5) 


N/A 


0.039 


0.337 


(0.01,0.04) 


(0.01.N/A) 


0.2420(6) 


N/A 


0.705(16) 


-0.162 


0.194 




(0.01,0.0343) 


N/A 


0.3328(6) 


N/A 


-0.035 


0.280 




(0.01,0.04) 


N/A 


0.3510(6) 


N/A 


-0.011 


0.295 


(0.02,0.04) 


(0.02.N/A) 


0.3228(6) 


N/A 


0.790(10) 


-0.218 


0.137 




(0.02,0.0343) 


N/A 


0.3681(6) 


N/A 


-0.142 


0.189 




(0.02,0.04) 


N/A 


0.3849(6) 


N/A 


-0.114 


0.208 


(0.03,0.04) 


(0.03.N/A) 


0.3880(7) 


N/A 


0.912(11) 


-0.391 


-0.020 




(0.03,0.0343) 


N/A 


0.4003(6) 


N/A 


-0.364 


-0.000 




(0.03,0.04) 


N/A 


0.4160(6) 


N/A 


-0.330 


0.025 
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404 



150 
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of q 2 obtained by the above fitting. 

The form factors in the physical kinematics are calculated from the extrapolation or in- 
terpolation with momentum and quark mass. For the physical kinematics of proton decay 
into meson and lepton final state, — q 2 is equivalent to lepton mass squared in the relevant 
form factor Wo(q 2 ). In the lattice computation, however, the quark masses are other pa- 
rameters that need to be tuned toward the physical pion and kaon masses. Therefore we 
have all together three parameters to tune: degenerate u, d quark mass m u d, strange quark 
mass m s and meson momentum \p\. In our simulation, the fh u d — > rh^/ s limit is taken by 
an extrapolation, fh s — > fh^ hys limit is taken by an interpolation, where physical quark mass 
is realized by the limit, 



m ud -> mfj s = 0.001385, 



iid 



m K 



_► mf ys = 0.03785, 



(31) 
(32) 
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FIG. 1: Effective mass plot of nucleon (top), pion (middle) and Kaon (bottom) at momen- 
tum square rip = (circle), rip = 1 (square), rip = 2 (diamond) which correspond to p = 
(0, 0, 0), (vr/12, 0, 0), (vr/12, 7r/12, 0) respectively. For nucleon we use gauge-invariant Gaussian 
source/sink, and for meson we use (Kuramashi-)wall source and gauge-invariant Gaussian sink. 
This is for the lightest quark mass m u d = 0.005 and m s = 0.0343. Solid line (colored band) 
indicate the central value (statistical error) obtained by fitting. 

with the values to reproduce the experimental hadron mass ratios, m n /mQ and mx/inn, the 



pion and kaon mass over the mass of fi~ 



We employ two different procedures for taking the above limit. One is to use the global 
fit with a function that depends on both quark mass and q 2 , and thus Wq at physical point 
is straightforwardly obtained. The other is to sequentially take both limits; first q 2 — > and 
then take the quark mass to the physical point. In this procedure Wq at physical point is 
obtained by the second limit. In the next section we will show numerical results with these 
procedures. 
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FIG. 2: Wq for p — > 7r° decay channel is plotted as a function of operator time (t in Eq. (I2ip ). 
The proton source is located at t = 5, and the 7r° sink is at t = 27. Different symbols show the 
two different proton interpolating fields, which correspond to (u T Cj^d)u (open) and {u T C^^^,d)u 
(filled). The horizontal solid line indicates the central value of constant fit to the both plateaus in 
the range 13 < t < 20 simultaneously. The shaded area indicates 1-sigma error bound. 

C. Extrapolation to physical kinematics with global fitting 

The global fitting to obtain the form factor in the physical kinematics uses 



F Wo( m ud,q 2 ) = A + A 1 m ud + A 2 q 2 , 
F§Afh ud , rh s , q 2 ) = B + Bpy + B 2 m s + B 3 q 2 



(33) 
(34) 
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FIG. 3: Wq for p — > tt° decay channel is plotted as a function of operator time. Symbols are same 
as in Figure [2J 



with free parameters Ai and Bi. F^f 1 is used for the pion or r\ finial state, F^ Q for the 
kaon final state. This procedure is the same as that employed in the previous study 16 ]. 
We use four different quark masses, two different strange quark masses and the two lowest 
non-zero spacial momenta, and therefore the total number of data points is eight for n and 
rj or sixteen for the kaon final states. The results obtained with the global fit using all 
the data are shown in the second column in Table [Til It turns out that the simple linear 
function as described in Eq.( l33l) and (134p is in good agreement with the lattice data for all 
channels, which is indicated by the reasonable x 2 /dof (< 1.4). The fit results F^(rh^/ S , g 2 ), 
F^ (m^/ S , mf ys , q 2 ) at the physical masses are shown in Figs. H] and [5j 
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FIG. 4: q 2 dependence of W^(q 2 ) at all quark masses in the lattice unit. We plot the results at 
f^ud = 0.005 (circle), m U( i = 0.01 (square), m u d = 0.02 (upper-triangle) and m uc i = 0.03 (down- 
triangle). In the figure for K 0,+ , results at m s = 0.0343 represent open symbol and filled symbol 
at m s = 0.04. In addition we also show the global fit function in the physical quark mass as solid 
line and its error band. The star symbol is a result in the physical kinematics using global fit. 
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FIG. 5: q 2 dependence of W^q 2 ) at all quark masses. Symbols are same as Figure HI 

D. Extrapolation to physical kinematics with sequential fitting 

In this procedure we first take the linear extrapolation or interpolation to q 2 = with two 
spatial momentum points in each mass m and second take a chiral extrapolation to physical 
quark mass. Figure |6] and [7] plot the results at q 2 = point as a function of fh u d after taking 
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the q 2 = limit. In the chiral extrapolation of the data at q 2 = the fitting function is 

fwo{™ud) = a + axrhud, (35) 

fw (™ud,rh s ) = b + b 1 m ud + b 2 m s , (36) 

for the pion, t] final state or kaon final state respectively. Here a^ and bi are the free fitting 
parameters. From Figure [6] and [7] we observe that the linear function well describes the 
lattice results for each matrix elements in four different mass points. Actually the x 2 /dof 
for all matrix elements are satisfied with x 2 /dof < 2. The results are shown in Table HT1 (see 
the column marked as "Sequential"). 

E. Systematic errors 

The systematic errors due to unphysical quark mass and kinematics used in the simulation 
(combined mass and q 2 limit), finite volume and non-zero lattice spacing will be discussed in 



this section. This work uses the lattice scale estimated in Ref. 18] and the renormalization 
constant shown in Eq. (f29]) and Eq. fl30l) . To yield the total error apart from the statistical 
error of the form factor we add in quadrature the systematic errors in the extrapolation, 
finite size effect and lattice artifact together with the error of lattice scale and of the non- 
perturbative renormalization procedure. 

At the target mass and momentum point (m ud ,m s ,q 2 ) = (m",mf ys , 0), no chiral 
singularity is expected. Therefore, if the simulations are made closer to the target, the 
linear approximation to the fitting function becomes arbitrary precise. However, as the 
simulated points might not be close enough to assume the linearity, we need to assess the 
systematic error due to the choice of this approximation. This systematic error is regarded 
as the effect of higher order than 0(rh u d) and 0(q 2 ). Note that the higher order effect than 
0(m s ) is safely neglected as its variation around the physical point is very small as can be 
estimated by comparing the results with m s = 0.0343 and 0.04 in Figs. |4]and|5j 

As the main results of the relevant form factors we employ those by the global fit with 
0.005 < m uc i < 0.03 (see in the second column of Table |TT])- For the convenience of later 
discussion we define the different fitting ranges as 

r full : [0.005, 0.03], r heavy : [0.01, 0.03], r light : [0.005, 0.02] (37) 
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FIG. 6: Results of Wq*(0) at different fh = m u d + m Ics . The different open symbols shown in 
the matrix element of Kaon final state are the results at different partially quenched strange quark 
mass m s = 0.0343 (circle), m s = 0.04 (square). Straight lines show linearly fit function with all 
four quark masses. For the matrix element of p — > K, these are the results after taking the physical 
strange quark mass. The cross symbol is result at physical light and strange mass with four fitting 
points and star symbol is with three fitting points. 
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FIG. 7: Results of Wq(0) with same symbols as Figj6j 

which are also used in Table [Til The variations of results removing the largest fh u d from the 
global fit, removing the smallest m u d from the global fit and the result in sequential fit from 
the main result provide the systematic errors coming from uncertainty of the fitting function 
for the extrapolation to the physical kinematics and finite size effect (FSE). 

The uncertainty in the extrapolation due to higher order effect than linearity of quark 
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mass (and also q 2 ) is estimated by variance between results in r fu11 and r llght and variance 
between results with global fit and sequential fit. By comparing the region with and without 
heavy mass m = 0.03 which is the strange quark mass, we estimate the 0(rh 2 ) effect. Fur- 
thermore since sequential fitting procedure, explained in the previous subsection, takes into 
account the mass-dependence of q 2 slope, we estimate the systematic error of the extrapo- 
lation to the physical kinematics as a part of the higher order effect, e.g. 0(rhq 2 ) terms, 
beyond the rh and q 2 linear approximation by comparing with results in the global fit. 

On the other hand the estimate of FSE using variance between results in r fu11 and r heavy 
expects to probe at least a part of FSE since the lightest point is affected most from the 
FSE rather than higher order mass effect. Such estimate of FSE has been known in the 



calculation of the nucleon axial charge qa 2% |25] in which significant FSE seemed to be 



observed in the lightest quark mass in the same gauge ensemble. (Furthermore this is also 
suggested noting that the relevant form factor Wo for a pion final state is proportional to 



1 -f- gA in the leading order of baryon chiral perturbation theory, see Ref. 16]). Therefore 



neglecting data at the lightest mass m = 0.005 from the fitting region might include less 



contamination of FSE (see also Fig. 10 of Ref. 251]) . 

As a consequence the systematic error including both 0(fh 2 ), 0((q 2 ) 2 ), 0(mq 2 ) and FSE 
is evaluated by the quadrature of variances of (global, sequential) fitting results in the range 
of r fu11 and maximum variance between global fitting results in the range of (r fu11 , r llght ) and 
(r , r heavy ) although it might be conservative. The magnitude of the above is shown in the 
column denoted as "Extrapolation" in Table IIHI 

The discretization error of 0(a) may arise from the inexact chiral symmetry due to finite 
L s , which is of order m rcs ~3x 10~ 3 and thus safely neglected. Here the dominant dis- 
cretization error at 0(a 2 ) has been estimated using the scaling study of hadronic observable 



performed with this and finer lattice ensembles [181 ] . The observed discrepancy in the spec- 
troscopy of light meson (Fig. 26 in Ref. 18[) with the two lattice spacings is up to 1-2 %, 
which amounts to about 5% discretization error for the form factor Wo assuming the 0(a 2 ) 
scaling. We take this 5% as the 0(a 2 ) error, whose magnitude appears as the naive estimate 
using the power counting (oAqcd) 2 ~ 0.02 with Aq CD = 250 MeV. 

In addition to all above, we take into account the error coming from uncertainty of lattice 
spacing which is given in error of a -1 = 1.73(3) GeV and the one for the renormalization 
constant which are given in Eq.f l29p and ( 13"U|) . 
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TABLE II: Table of results for renormalized W (fi = 2GeV) in GeV 2 after global and sequential 
fitting. The error is only statistical one. For global fitting, we show the results with three different 
fitting mass-ranges, which are all range 0.005 < m uc i < 0.03 (r fu11 ), exclude the heaviest mass, 
m uc i = 0.03, (r llght ) and exclude the lightest mass, m u d = 0.005, (r heavy ). For the sequential fitting, 
we show the results with all mass range. 







Global 




Sequential 


matrix element 


r full 


*7dof 


flight 


~, heavy 


„,full 


X 2 /dof 


(7T°\(ud) R U L \p) 


-0.103(23) 


1.4 


-0.132(29) 


-0.072(34) 


-0.114(22) 


2.2 


(7T°\(ud) L U L \ P ) 


0.133(29) 


1.4 


0.156(41) 


0.142(38) 


0.123(28) 


1.1 


(K°\(us) R u L \p) 


0.098(15) 


0.4 


0.103(19) 


0.092(29) 


0.093(15) 


0.1 


(K°\(us) L u L \p) 


0.042(13) 


0.4 


0.044(16) 


0.037(20) 


0.037(14) 


0.1 


(K+\(us) R d L \p) 


-0.054(11) 


0.8 


-0.060(13) 


-0.052(21) 


-0.049(13) 


0.6 


(K+\(us) L d L \p) 


0.036(12) 


0.8 


0.040(15) 


0.041(18) 


0.041(12) 


0.6 


(K+\(ud) R s L \p) 


-0.093(24) 


0.6 


-0.108(28) 


-0.082(39) 


-0.088(25) 


0.9 


(K + \(ud) L s L \p) 


0.111(22) 


0.6 


0.121(28) 


0.115(37) 


0.117(23) 


0.7 


(K+\(ds) R u L \p) 


-0.044(12) 


0.1 


-0.043(14) 


-0.041(20) 


-0.044(12) 


0.1 


(K+\(ds) L u L \p) 


-0.076(14) 


0.3 


-0.082(17) 


-0.076(24) 


-0.078(14) 


0.5 


(r]\(ud) R u L \p) 


0.015(14) 


1.3 


-0.002(19) 


0.031(19) 


0.017(14) 


1.2 


(r)\(ud) L u L \p) 


0.088(21) 


0.7 


0.094(29) 


0.094(28) 


0.076(21) 


0.4 



We ignore the partially quenched effect of strange quark, which is due to the small 
mismatch of the sea and valence strange masses, for the matrix element of K meson final 
state. Since the valence strange quark mass dependence of Wq is negligibly small as shown 
in Figj6]and FigJTJ this effect is also negligible. Note that we also do not consider the effect 
of disconnected diagrams in the matrix elements of the r\ in the final state, but note that 
the result is valid assuming flavor SU(3) degenerate valence quark vn^\ = mj al and ignoring 
partially quenched effect of strange quark. 
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F. Results of proton decay matrix elements 

Table ITTT1 summarizes the results of the relevant form factor Wo(q 2 ) of proton decay for 
all the principal matrix elements Eqs. ( ITS"]) . ( ir5|) - (fTU|) at q 2 = 0. The central values are 
those obtained with the global fit on q 2 and the simulated quark masses for the physical 
kinematics rh u d — > fn^d* ■> "fhs —*■ rh^ hys and q 2 — ¥ 0, with the r fu11 range for m U d, which has 
already been shown in Table [TT1 Values in the first parentheses are the statistical errors. The 
systematic error budget is shown in the last four columns. These four errors are added in 
quadrature to give the total systematic error shown in the second parenthesis for each form 
factor value. 

Figure [H] shows the results of the form factors with the error bars expressing the total error 
of statistical and systematic errors added quadrature, which are marked as "iV/ = 2 + 1". 
The two panels compare the results with old ones using some approximation. The left panel 
compares against the results with quenched approximation in direct method [16|. The right 



panel shows those with the indirect method in the same ensembles 17|. The sizable error for 
u Nf = 2 + 1" in the current analysis prevents us from seeing any visible difference from the 
quenched or indirect results. For phenomenological application, however, one should clearly 
use the Nf = 2 + 1 results with their total error instead of the previous results, because each 
approximations lead to an unspecified systematic uncertainty. 

V. SUMMARY AND OUTLOOK 

We have presented the lattice calculation of proton decay matrix elements using 2 + 1 
flavor dynamical domain-wall fermions, essential ingredients to reliably estimate the nucleon 
lifetime in grand unified theories. The direct method using three-point function (nucleon)- 
( operator)- (meson), with non-perturbative renormalization, was applied on a volume L\ ~ 3 
fm 3 . Previous calculations had unestimated systematic uncertainties on the matrix elements 
at the physical kinematics. This work made it possible to control these uncertainties for the 
first time, by removing most of them, while remaining uncertainties were given with their 



estimates. The uncertainties t 



Q 



rat have been eliminated are those due to the quenched ap- 



proximation [16| and the use 17] of the indirect method with the tree-level baryon chiral 



perturbation theory. The estimated uncertainties are the combined errors from chiral extrap- 
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J I R 

TABLE III: Final results of renormalized W (fi = 2GeV) for individual matrix elements and 

T / R 

error budget of statistical and systematic uncertainties. The first and second errors in W 
represent statistical and systematic ones respectively. The third column denotes total error which 
is estimated by adding in quadrature statistical and systematical errors. The fourth column denoted 
as x shows the systematic error of mass and momentum extrapolation/interpolation estimated by 
the variance of extrapolation to physical kinematics and fifth column is uncertainties from lattice 
artifacts explained in the text. The last two columns show the uncertainties of renormalization 
factor (AZ) and lattice spacing (Aa _1 ). We also show the p — > 7r + u decay matrix element using 
SU(2) isospin relation in Eq. ([20|) . 











Total error Systematic 


3rror budget 


Matrix element Wq(h = 


2GeV) GeV 2 


(%) 


X 


0(a 2 ) 


AZ Aa" 1 


(7T°\(ud) R UL\p) 


-0.103 


(23) 


(34) 


40 


0.033 


0.005 


0.008 0.004 


(TT°\(ud) LUL \p) 


0.133 


(29) 


(28) 


30 


0.026 


0.007 


0.011 0.005 


(7T+\(ud) R d L \ P ) 


-0.146 


(33) 


(48) 


40 


0.047 


0.007 


0.011 0.006 


(7T + \(ud) L d L \p) 


0.188 


(41) 


(40) 


30 


0.037 


0.010 


0.016 0.007 


(K°\(us) R u L \p) 


0.098 


(15) 


(12) 


20 


0.007 


0.005 


0.008 0.003 


(K \(us) L u L \p) 


0.042 


(13) 


.8) 


36 


0.007 


0.002 


0.003 0.001 


(K+\(us) R d L \p) 


-0.054 


(11) 


.9) 


26 


0.008 


0.003 


0.004 0.002 


(K + \(us) L d L \p) 


0.036 


(12) 


(7) 


39 


0.007 


0.002 


0.003 0.001 


(K+\(ud) R s L \p) 


-0.093 


(24) 


(18) 


32 


0.016 


0.005 


0.008 0.003 


(K + \(ud) L s L \p) 


0.111 


(22) 


(16) 


25 


0.012 


0.006 


0.009 0.004 


(K+\(ds) R u L \p) 


-0.044 


(12) 


(5) 


30 


0.003 


0.002 


0.004 0.002 


(K + \(ds) L u L \p) 


-0.076 


(14) 


0) 


22 


0.006 


0.004 


0.006 0.003 


(v\(ud) R u L \p) 


0.015 


(14) 


(17) 


147 


0.017 


0.001 


0.001 0.001 


(r]\(ud) L u L \p) 


0.088 


(21) 


(16) 


30 


0.014 


0.004 


0.007 0.003 



olation and finite volume, discretization error, error in the non-perturbative renormalization 
and the estimate of the lattice scale. The relevant form factors Wo(q 2 = 0) of the twelve 
principal matrix elements Eqs. (TT3|) . ( TT5|) -( TT9|) . from which one can calculate those for all 
the nucleon to pseudoscalar meson process, has been evaluated and summarized in Table III II 
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FIG. 8: Summary of W Q (/i = 2GeV) for twelve principal matrix elements. Filled circles show 
the present results, and for the comparison the results in quenched QCD (open circle) and indirect 
method using chiral perturbation theory (cross) are plotted in the same raw. 

with their error estimates. 

In this paper we have established an estimate of proton decay matrix element with all the 
errors and note that the errors are fairly large (30%-40 % for ir final state and 20 %-40% total 
error for K final state). One of the major uncertainty is the statistical error, especially for 
p —?• e + 7r° decay mode, and that could have influenced the size of the error of combined chiral 
extrapolation and finite volume effect. A significant improvement of the current results is 
expected by adopting the newly developed technique for reduction of the statistical error 



26| . which will be addressed in future work. We want to emphasize, though, for now in 



any serious phenomenological application one should use the results in this study with their 
total errors. 
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